Phonons in Random Elastic Media and the Boson Peak 
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We show that the density of states of random wave equations, normahzed by the square of the 
frequency, has a peak — sometimes narrow and sometimes broad — in the range of wave vectors 
between the disorder correlation length and the interatomic spacing. The results of this letter may 
be relevant for understanding vibrational spectra and light propagation in disordered solids. 
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One of the intriguing features of random elastic media 
— observable in both Raman or neutron cross sections 
but also in calorimetric measurements — is an anoma- 
lous accumulation of phonons at low frequencies 0] . This 
phenomenon finds its most prominent manifestation in a 
peak in the quantity = p{uj)/llP', where p{uj) is the 
phonon density of states (DoS) and lo is the frequency. 
In recent years, many competing theories as to the origin 
of this so-called boson peak have been formulated (for a 
recent reference, see Ref. and references therein). Of 
these approaches a majority is based on model mecha- 
nisms specific to the low temperature physics of amor- 
phous solids. Yet the boson peak is shown by both 
glassy and random crystalline systems alike, an obser- 
vation which has ignited the search for an explanation 
which is not tied to the specifics of a glassy environment. 



On a basic level, the acoustic excitations of both amor- 
phous materials and disordered crystals are described by 
random wave equations. Existing analyses of such equa- 
tions in the literature indeed predicted a disorder gen- 
erated excess DoS, 0. However, these structures where 
observed at high frequencies (wavelengths of the order of 
the interatomic spacing), while the boson peak is a low 
energy phenomenon. In this letter we argue that the DoS 
of elastic vibrations in disordered media is enhanced by 
disorder at low frequencies, with /(w) exhibiting either a 
peak or a broad maximum (see Fig.nfor a- representative 
picture of the DoS). Which of the two will be observed 
crucially depends on the modelling of the randomness, i. 
e. the type of disorder at work. Since the profile of the 
randomness of 'real' systems is generally unknown, we 
are not in a position to judge whether the present mech- 
anism alone may account for the spectral peaks observed 
in all experiments. We believe, however, that it operates 
under quite general conditions. 



Let us begin by discussing some generic large scale 
structures of the DoS of acoustic excitations. In a clean 
system, the linear relation u> = ck between frequency and 




FIG. 1: for the random elastic media in cases where 

0:2 > or 02 < 0. 



magnitude k of the wave vector gives rise to a DoS 
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Here, Ad is the area of the d-dimensional unit sphere 
and c the speed of sound (for our present analysis, the 
potential existence of several acoustic phonon branches 
with different sound velocities will not be of importance). 
Tied to the linearizability of the dispersion relation, the 
relation p uj'^~^ applies only to frequencies oj ^ uju 
much smaller than the Debye frequency cod ~ f (a: lat- 
tice spacing.) On dimensional grounds, the generaliza- 
tion of 1^ to higher frequencies must be of the form 
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where f(u) is some function with f{u <C 1) ~ 1. It is 
clear that f{u) falls off to zero at u ^ 1. 

In contrast, in disordered materials we find strong de- 
viations from the above scaling form already for frequen- 
cies CO <ti LU£). Indeed, our analysis of the random wave 
equation below will lead us to 
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where the effective correlation length ^ ^ a of the disor- 
der is assumed to be much larger than the lattice spacing, 
and c is the typical speed of sound in random media to 
be defined more precisely below. The scaling function / 
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is given by 

f{u) = 1 + aiu'^ + u < 1, 

f{u) = + ^ + ...), (4) 

where the form of the coefficients ai^2, P depends on the 
space dimensionality and on a few basic characteristics 
of the randomness. In the most interesting case, d = 3, 
we find that /3 > 1 and ai > 0. The sign of a2 lacks 
universality. a2 is negative if disorder in materials is 
mostly in the random elastic constants, and it is positive 
if disorder is mostly due to fluctuating mass density. 

For a2 > 0, grows for u <C c/l and falls off at 

uj ^ c/l (while still lo <ti ldd). This implies the existence 
a peak at w ~ c/l. In contrast, for a2 < 0, increases 
until w becomes much larger than c/l. Combined with 
the expected drop-off of /(w) at w > lojj, this produces 
a broad maximum for I{lu) between c/l and ujd- Fig- H 
shows a caricature of the two scenarios. We finally note 
that if the disorder is relatively weak, rms c < c, (so that 
the fluctuations of local speed of sound c are smaller than 
the typical sound velocity) ai ~ a2 ~ var c/c^; for more 
generally applicable expressions, see below. 

We next turn to the derivation of these results. Con- 
sider the wave equation 

where the random velocity field c(x) fluctuates on spatial 
scales ^ I. Interpreting the variable to(x) = l/c^(x) 
as the mass density of the random medium, we refer to 
Eq. ISJ as a 'random mass' wave equation. 

In the regime of small frequencies, lo <C c/l, we 
are facing a situation where the correlation length I 
is much smaller than the 'typical' wave length = 
c/lu at which the wave function V'(x) fluctuates. The 
phonon field effectively averages over many fluctuation 
intervals of the disorder and, to a first approximation, 
Eq. lO may be replaced for its average over random c(x): 

+ ipi'^) ~ 0, where the typical speed of sound c 
is defined as 

c^(c(x)-2)-V2 (6) 

and the angular brackets denote averaging over random 
c(x). The very low frequency asymptotics is then given 
by Eq. |^ with c taken from Eq. ©. (This result, as 
well as the high frequency asymptotics Eq. (|f 0|l discussed 
below, are originally due to Chalker Q.) 

To compute corrections to the low frequency asymp- 
totics of the DoS we rewrite Eq. (0 as 

-A-^^Mx))^(x) = ^^^(x), (7) 

where the function /i(x) = [c/c(x)]^ — 1 describes the ran- 
domness. This representation suggests to interpret uP' /l? 



as the eigenvalue of the operator A weakly perturbed by 
Lo^h{-x)/<?. The unperturbed problem is trivially diago- 
nalized by a set of plane waves V'k(x) = exp(ik • x), 
with eigenvalues uP' /<? = k"^ {L is the system size.) We 
next apply standard perturbation theory to compute the 
eigenvalue shift caused by the presence of the pertur- 
bation uPh/(?. The vanishing of (ft.) implies that, on 
average, there are no first order corrections. To sec- 
ond order in h we find that the average eigenvalue is 
given hy ^ ^e + k^j (Cf^ fc^-fc'^ ^ ^^ere ^(k) is the 
Fourier transform of the disorder correlation function of 
(/(x — y) = ( h{x.)h{y) ). Irrespective of the distribution 
of the disorder, this function (a) drops off rapidly for 
k > l^^ and (b) approaches a constant value ~ l'^var{h) 
for k ^ l^^. This implies that for d > 2, the integral 
above is dominated by momenta k' l^^. At the same 
time, the reference momentum k ~ uj /c <^ l~^ is small. 
We thus neglect the ^-dependence of the integrand and 
arrive at 

Physically, the correction to the zeroth order eigenvalue 
is due to virtual scattering events wherein states with 
low-lying momentum k lu/c scatter off rapid fluctua- 
tions of h into high-lying states k' ~ l~^. Yet for small 
frequencies, the large phase volume ~ l~'^ available to 
these scattering processes cannot outweigh the overall 
multiplicative factor oj^. This mechanism pervades to 
higher orders in the eigenvalue expansion and justifies 
the perturbative approach. In particular, Eq. ^ indeed 
describes the dominant correction to the low frequency 
dispersion relation. 

We next substitute Eq. ((HJ into p{uj) — 
-j^j (5(cL;(fc) — Lu) to obtain an expansion of the DoS 
as in Eqs. (|3I4|) . Specifically, the coefficient 

d + 2 f d'^k g(k) 
-i = ^y (^%^-var(ft)>0, (9) 

where the proportionality to var(ft,) follows from the 
properties of the correlation function g discussed above. 
Interestingly, positivity of ai in Eq. Q is a direct con- 
sequence of the fact that the second order perturbation 
theory always lowers the ground state energy. 

For small velocity fluctuations, var(c) <C c^, we obtain 
the estimate ai ~ var(/i) ~ var(c)/c^ quoted above. We 
finally note that these results hold only for d > 2. For 
d < 2, the integral in Eq. Q is infrared divergent and 
the approximation scheme employed here breaks down. 

We next turn to the discussion of the high frequency 
case, (jj ^ c/l. In this regime, the velocity field varies 
very little over length scales comparable to the typical 
wavelength. This implies that, locally, the solutions of Q 
behave like plane waves with the local dispersion relation 
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fc(x) = (jj/c(x) and the local density of states given by 

Ad 



UJ 



d-1 



Averaging this result over the random velocity field we 
obtain the global density of states 
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We thus find that the coefficient /?, introduced in Eq. Q), 
is given by 



(C-2(X))^ 



(10) 



As a consequence of the convexity of the power function, 
/3 > 1 for d > 2. 

To obtain corrections to the high frequency asymp- 
totics (jlO|l . we need to compute distortions in the spectral 
density caused by shallow (compared to the wave length) 
variations of the velocity field. It turns out that this task 
is most efficiently tackled by analyzing the Green's func- 
tion 



m(x) 



(11) 



where m = c ^. From ()11|) the average DoS is obtained 

as 



Im(m(x) G(x,x)) 



(12) 



What makes the operator pi|l a good starting point for 
our analysis is its structural similarity to a Schrodinger 
operator with 'Planck's constant' h ^ oj^^. This anal- 
ogy will enable us to apply semiclassical approximation 
schemes familiar from quantum mechanics. We begin by 
applying the Wigner transform 

G(x,y) = ^ J d^k G(x,k) e*-k.(x-y)^ 
whereupon the 'Schrodinger equation' assumes the form 



fk — iuj 



m(x) 



G(x,k) = 1. 



(13) 



To make use of the smallness, ^ oj ^, of the derivative 
operators, we expand G(x, k) in powers of lo~^ , 

G(x, k) = G(") (x, k) + -G^^^ (x, k) + ^^G^^) (x, k) -f . . . , 

UJ UJ'^ 

substitute the expansion into Eq. H13|) . and find G*^°-', 
G'"'^', G*^^-', . . ., recursively. By symmetry, G^^-' vanishes 
upon configurational averaging so that the dominant cor- 
rection to the DoS is provided by G''^-'. Substituting this 



coefficient into Eq. 112() , and comparing with Eq. we 
obtain the result 



a2 = 
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Clearly, a2 > if = 3. If the fluctuations of c are weak, 
a2 ^ var c/c^ ~ ai. 

We thus derived Eq. Q). The positivity of the two ex- 
pansion coefficients ai_2 implies that the scaling function 
f{u) in Eq. JSJ grows at u <C 1 and falls off at u ^ 1. 
This implies that the functional profile of the DoS con- 
tains a peak somewhere at u = 0(1). 

To conclude our analysis of the prototypical wave equa- 
tion (0) , let us briefly discuss the case of space dimensions 
d 7^ 3. For d — 1^ the high frequency coefficient q;2 > 0, 
yet /3 < 1. Due to the instability in the the low frequency 
expansion, the methods applied in this letter do not al- 
low to calculate ai at d < 2. However, for certain dis- 
tributions of the disorder transfer matrix methods may 
be applied to obtain an exact solution 0]. These calcu- 
lations show that f{u) is a globally decreasing function, 
and ai < 0. For d ~ 2, the two dominant high-frequency 
coefficients are structureless, /3 = 1 and 012 = 0, while the 
low frequency expansion in Eq. Q will now involve terms 
proportional to log(u). From these results we cannot de- 
cide whether p(uj)/uj has a maximum or a minimum be- 
tween its low and high frequency asymptotics. Finally, 
for d > 3, Q!i > and (3 > 1, however a2 < 0, which 
implies that f{u) is a monotonously increasing function. 
Summarizing, we see that the normalized DoS of Eq. jsj, 
p{uj)/uj'^~^, exhibits a peak at wavelengths of the order 
of disorder correlation length only for d — 3. 

Before proceeding, let us briefly compare our so far re- 
sults to earlier work. In most numerical simulations of 
the problem (cf. e.g., Ref. and references therein), the 
disorder is chosen to be uncorrelated, that is, its corre- 
lation length is of the order of the lattice spacing. For 
such type of disorder, m ~ 1 translates to uj ud deep 
in the bulk of the spectrum. At these frequencies it is 
hard to tell whether deviations from the low frequency 



asymptotics p{uj) 



are caused by lattice effects (cf. 



Eq. ^) or by disorder (Eq. (jSJ). Qualitatively, however, 
the numerical data is in agreement with the results of our 
present analysis. Turning to earlier analytical work, we 
notice that most approaches to random elastic problems 
rely on the self consistent Born approximation (SCBA). 
(See Ref. |G||for a general review of the methods involved 
and Refs. [3, 01 for the appfications of these methods 
to random wave equations.) However, for a number of 
reasons the results obtained by this method lack quan- 
titative reliability when used to calculate the DoS: (i) 
Within the standard SCBA, 1/c^ is modelled as a Gaus- 
sian distributed variable. This implies the existence of 
rare domains where 1/c^ is negative. Using the tech- 
niques introduced in Ref. Q it is possible to show that 
these domains, no matter how small, lead to unstable 
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wave modes and to unphysical sharply growing contri- 
butions to the density of states at higher frequencies, 
(ii) The SCBA neglects certain contributions ('crossed 
diagrams') to the perturbative expansion of the Green 
function At high frequencies, this approximation 

becomes invalid. The value of the threshold frequency 
beyond which the SCBA breaks down depends on the 
type of disorder under consideration. At any rate, how- 
ever, it is parametrically smaller than the band width. At 
the same time, (iii) , previous studies focused on the case 
of short range correlated disorder. In this case, deviations 
of the scaling law p ~ lo'^^^ are expected at frequencies 
Lo uj£), i.e. in a regime where the SCBA breaks down 
and lattice effects intervene. (The alternative techniques 
in this letter avoid all of these problems.) 

The wave equation Eq. (jSJ applies to the specific case, 
where only the mass density of the elastic medium fluc- 
tuates. The continuum description of a more general 
environment wherein the elastic constants also fluctuate 
reads as 



(V/i(x)V + cj2to(x)) V'(x) 0. 



(14) 



Here both /j.(x) and m(x) are random positive quantities. 

It turns out that the direct perturbative expansion ap- 
plied above to Eq. © cannot be used to determine the 
coefficient ai of the problem Eq. H14|l . As an alternative, 
we apply a generalized variant of the self-consistent Born 
approximation, wherein /i(x) — ^jLq -\- (T^(x), and cr(x) is 
a Gaussian distributed variable with zero mean. (In this 
way, positivity of the elastic constant is ensured.) The 
actual implementation of the SCBA for this type of dis- 
order turns out to be rather involved and its details will 
be published elsewhere. Suffice to say that at d = 3, 
ai is still positive. a2, on the other hand, no longer 
has definite sign. Applying the high frequency expansion 
outlined above for the random mass density case we find 
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If Q!2 is negative (which happens, for example, in the 
limiting case of non-random to(x)) f{u) no longer has 
a maximum at u ~ 1. Instead, I{uj) is an increasing 
function of lo with a broad maximum reached at frequen- 
cies u > c/l, before dropping off at frequencies higher 
than ujd. Concluding, we find that, depending on the 
profile of the disorder, the envelope function of the 
random mass/elastic constants problem either contains 
a low-frequency peak, or a broad high frequency maxi- 
mum. 

So far, we have considered the case of scalar phonons. 
In a realistic environment, however, ip —>■ Ui will be a 



d-component vector. The most general random vector 
problem would be governed by a formidable rank four 
random elastic modulus tensor. Assuming, however, a 
medium consisting of a random accumulation of 'micro- 
crystallites' each of which possessing intrinsic rotational 
invariance — this assumption may well be violated, espe- 
cially in glassy environments, although it seems to work 
in polycrystalline materials — the effective wave equation 
reduces to 

di [X{x)djUj] + dj [/i(x) {diUj + djUi)] + Lj'^m{x)ui = 0, 

(16) 

where m(x) is a random density of the medium and A(x), 
/i(x) are random Lame coefficients. 

In the limit where only A(x) is random, Eq. H16() can be 
mapped into Eq. Q by substitution of ?A = diUi{X + 2p). 
In this case, both ai and a2 are positive. This result 
carries over to the case where, in addition to random 
A(x), also the density to(x) is made random; For these 
two types of disorder I{uj) has a peak at w ~ c/l. If, 
however, the shear modulus p is also random, Eq. H16(l 
becomes more similar to Eq. (|14f) . Applying the SCBA, 
it is still possible to show that ai > (for d > 2). As 
with our previous discussion of H14|l . however, a2 is no 
longer of definite sign. Specifically, I{uj) will contain a 
broad maximum (q!2 < 0) if only the shear modulus p{x) 
is random. 

To conclude, we have shown that, depending on the 
type of disorder, the normalized DoS I{uj) of random 
wave equations may either contain a peak at phonon 
wavelengths of the order of the disorder correlation 
length, or a broad maximum at wavelengths below the 
correlation length. We believe that this work is not 
only relevant for the interpretation of data on vibrational 
modes in random media but also to the analysis of other 
types of waves, such as electromagnetic waves in random 
environments. 
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